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Collective Dynamics in Arrays of 
Coupled Nonlinear Resonators 

Ron Lifshitz, Eyal Kenig, and M.C. Cross 



11.1 Arrays of Nonlinear MEMS & NEMS Resonators 

The study of collective nonlinear dynamics of coupled mechanical resonators is re- 
gaining attention in recent years thanks to rapid developments in the fields of micro- 



electromechanical and nanoelectromechanical systems (MEMS and NEMS) (Roukes, 



2001[|Cleland, 2003[ ). MEMS & NEMS resonators are typically characterized by very 
high frequencies, extremely small masses, and weak damping. As such, they are natu- 
rally being developed for a variety of applications such as sensing with unprecedented 



accuracy ( 


Rugar et al, 2004 


Ilic et al, 2004 


1 


Ekinci et al, 2004 


Yang et al, 2006 


Li et al, 2007 


Naik et al, 2009 


Lee et al, 201C 


)). NEMS, in particular, are being de- 



veloped also for studying fundamental physics at small scales — exploring mesoscopic 



phenomena ( 


Schwab et al, 2000 


Weig et al, 2004 


I, and even approaching quantum be- 


havior ( 


LaHaye et al, 2004 Naik et al, 2006 


Rocheleau et al, 2010 


O'Connell et al. 



2010). MEMS & NEMS resonators often exhibit nonlinear behavior in their dynamics. 



as recently reviewed by Lifshitz and Cross (2008 2010) and by Rhoads et al (2010) 



This includes nonlinear resonant response showing frequency pulling, multistability, 



and hysteresis ^Craighead, 2000 Turner et al, 1998 Zaitsev et al, 2005 Aldridge and 



Cleland, 2005 Kozinsky et al, 2007), the appearance of chaotic dynamics (Scheible 



well as the formation of extended (Buks and Roukes, 2002) and localized (Sato et al. 



2006 Sato and Sievers, 2007 Sato and Sievers, 20081 collective states in arrays of 



et al, 2002[ |DeMartini et al, 2007"[ [Karabalin et al, 2009} [Kenig et al 2011[), as 



coupled nonlinear resonators. Nonlinearities may be a nuisance in actual applications, 
and schemes are being developed to avoid them, as demonstrated, for example, by 



Kacem et al (2009 20101. On the other hand, one can also benefit from the existence 



of nonlinearity, for example in mass-sensing applications (Zhang et al, 2002 Buks and 



Yurke, 2006), in suppressing noise induced phase diffusion, as suggested by Greywall 
et al (1994), and in achieving self-synchronization of large arrays, as proposed by 



Cross et al ( |2004[ |2006[ ) . Nonlinearity is even proposed by Katz et al (|2007| [2008| as 
a way to detect quantum behavior in large mechanical systems. 

Current technology enables the fabrication of large arrays, composed of hundreds 
to tens of thousands of MEMS and NEMS devices, coupled by electric, magnetic, or 
elastic forces. These arrays offer new possibilities for quantitative studies of nonlinear 
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dynamics in systems with an intermediate number of degrees of freedom — much larger 
than one can deal with in macroscopic experiments, yet much smaller than one con- 
fronts when considering nonlinear aspects of phonon dynamics in a crystal. Our studies 
of collective nonlinear dynamics of MEMS and NEMS were originally motivated by 



the experiment of Buks and Roukes (2002). These studies have led to a quantitative 
understanding of the collective response of arrays of nonlinear resonators, providing 
explicit bifurcation diagrams that explain the transitions between different extended 
modes of an array as the strength and frequency of the external drive are varied qua- 



sistatically (Lifshitz and Cross, 2003 Bromberg et al., 2006 1. We have considered more 



general issues such as the nonlinear competition between extended modes, or patterns, 
of the system — when many such patterns are simultaneously stable — as the external 



driving parameters are changed abruptly or ramped as a function of time (Kenig et al., 



2009a). We have also studied the formation, stability, and rich dynamics of intrinsi- 
cally localized modes (Kenig et al., 20096). Furthermore, we have investigated the 



synchronization that may occur in coupled arrays of non-identical nonlinear oscilla- 
tors, based on the ability of nonlinear oscillators to tune their frequency by changing 



their oscillation amphtude (Cross et al., 2004 Cross et al., 2006) 



The purpose of this chapter is to provide a review of the collective dynamical phe- 
nomena observed in these different systems, while highlighting the common concepts 
and theoretical tools that we have developed for dealing with them. We assume that 
the reader is familiar with the basic dynamical phenomena associated with single non- 
linear resonators. The unfamiliar reader is encouraged to consult our previous review 



2010) 



on the subjec t ([Lifshitz and Cross, 2008[ ), or its revised version ( [Lifshitz and Cross, 
In Sec 



11.2 we describe the equations of motion that are used to model arrays 



of nonlinear MEMS and NEMS resonators, for different experimental realizations. We 
then give two examples of the derivation and then application of discrete amplitude 
equations for treating arrays of resonators — in Sec. [11.3[ we study the resonant nonlin- 
ear response of arrays to parametric excitation, and in Sec. [11.4[ we discuss the question 
of synchronization. We conclude with two examples of the derivation and then appli- 
cation of continuous amplitude equations for treating large arrays of resonators — in 
Sec. [11.5[ for investigating pattern selection, and in Sec. [11.6[ for the study of intrinsi- 
cally localized modes. In place of a formal concluding section, we wish to emphasize at 
the outset that all the results obtained from analyzing the different amplitude equa- 
tions are in excellent agreement with numerical solutions of the underlying equations 
of motion. This upholds the validity of using such reduced descriptions for complex 
systems, whose original description is given in terms of coupled nonlinear ordinary dif- 
ferential equations. Furthermore, our numerical simulations of the equations of motion 
suggest that the predicted effects can be observed in arrays of real MEMS and NEMS 
resonators, thus motivating new experiments in these systems. 



11.2 Equations of Motion and Basic Assumptions 

Typical MEMS and NEMS resonators are characterized by extremely high frequen- 



cies — now going beyond 10 GHz (Huang et al, 2003 Cleland and Geller, 2004 Wein- 



stein and Bhave, 2010 ) — and relatively weak dissipation, with quality factors Q in the 
range of 10^ — 10^. For such devices, under external driving conditions, transients die 
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out rapidly, making it is easy to acquire sufficient data to characterize the steady-state 
well. This, and the fact that weak dissipation and weak nonlinearity can be treated 
perturbatively, are a great advantage for quantitative comparison between theory and 
experiment. 

11.2.1 Modeling a single nonlinear resonator 

A typical single resonator is described after appropriate scaling by a dimensionless 



equation of motion of the form ( Lifshitz and Cross, 2008 1 

X + Q^^i + [1 + H cosujpt] X + + -qx^x = G cos {uJot + <j)g) ■ (H-l) 

We typically use the fact that damping Q^^ is much smaller than the resonant fre- 
quency, which has been scaled here to 1, to define the small expansion parameter 
e — Q^^. The term proportional to H on the left hand side is an external drive that 
modulates the spring constant. This is a parametric drive — a term that is proportional 
to the displacement x as well as to the strength of the drive. The term proportional to 
G on the right-hand side is the standard direct drive, possibly shifted by a phase 0g 
with respect to the parametric drive. The coefficient of the nonlinear x^ Duffing term 
has been scaled to 1, and a nonlinear damping term [Dykman and Krivoglaz ( jl975| 



1984 )] with coefficient r] is also added. For parametric drive, we normally consider the 
largest excitation effect that occurs when the pump frequency ujp is close to twice the 
resonant frequency of the resonator. We therefore take Wp — 2 + eil.p, and take the 
drive amplitude to scale as the damping by setting H — eh. The amplitude of the 
direct drive is scaled as G = e^^^5, and its frequency is set an amount cQd away from 
the resonant frequency. 

The scaled equation of motion that we then obtain is of the form 

x + ex+{l + ehcos [(2 -I- enp)t]) x + x^ + f]x'^x = e^/'^\g\ cos [(1 + ef^n)t + (fig] , (11.2) 

where we use g = \g\e^'^!> to denote a complex drive amplitude. Ignoring transients, 
which as explained above decay very rapidly, the solutions to such equations of motion 
are of the form x — y/e 3?{v4(T)e**} plus corrections of higher order in e, where secular 
perturbation theory is used to yield the equation that governs the slow dynamics of the 
complex amplitude A(T). The variable T — et is the slow time scale upon which the 
interesting nonlinear dynamics takes place. Please refer to Lifshitz and Cross (20081 



2010 1 for more details and for many examples of the use of this approach. We only wish 



to remind the reader that additional nonlinear terms, up to third order in x or i such 



as x'^ and xx'^, that seem to be missing in eqn (11.2), merely conspire to renormalize 
the effective parameters in the slow equation for A{T), but do not affect the actual 
form of this equation. We therefore ignore all such terms as they have no effect on the 
actual nature of the solutions that we study. 

11.2.2 Modeling an array of nonlinear resonators 



Lifshitz and Cross (2003 1 had originally modeled a 1-dimensional array of parametri- 



cally driven coupled nonlinear resonators, motivated by the particular experiment of 



Buks and Roukes (2002), in which an array of 67 doubly-clamped micromechanical 
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gold beams was parametrically excited by modulating the strength of an externally- 
controlled electrostatic coupling between neighboring beams. We used a set of coupled 
equations of motion of the form 

Un +Un - ^Q^^{Un+l - 2u„ + U„_i) + ^ (D + HcOSLOpt) (?i„+i - 2u„ + U„_i) 

+ul ~ ^ri[{Un+l " W„)2(u„+i - ?i„) - (u„ - u„_i)^(u„ - iln-l)] = 0, (11.3) 

where u„(t) describes the deviation of the n*'' resonator from its equilibrium, with 
n = 1 . . . N, and fixed boundary conditions uq = un+i = 0. Detailed arguments for 
the choice of terms introduced into these particular equations of motion are discussed 



by Lifshitz and Cross (2003). We only note that they contain nearest-neighbor linear 
coupling which is both reactive, proportional to the relative displacements, and dissi- 
pative, proportional to the relative velocities; as well as nonlinear dissipative coupling, 
proportional to the square of the relative displacements and to the relative velocities. 
A simpler model, suitable in many other situations, is to take the equation of 



motion of each resonator to be as in ( 11.1 ) with the addition of only a linear reactive 



coupling term to its two neighbors. The equations of motion then take the form 

Un + Q^^Un + (l +H COS UJpt)Un+U^+T]u1un + ID{u„+i - 2m„ = 0, (11.4) 

where in both cases one could add the direct drive, proportional to G, that was con- 



sidered earlier in eqn (11.1) 



Finally, to model an array of oscillators — having a frequency-independent source 
of energy that sustains their oscillations — rather than simple resonators that respond 
resonantly to an external frequency-dependent drive, we consider a slight modification 



of eqn (11.4), given by 



it„ -I- Lolun - - U^n)Un + au^ + ^D{Un+l - 2li„ + M„_i) = 0. (11-5) 

In this case both the parametric drive and the direct drive are omitted. Instead, 
we introduce a negative linear damping with coefficient v, which represents an energy 
source to sustain the oscillations, while keeping the positive nonlinear damping so that 
the oscillation amplitude saturates at a finite value. We use a different scaling than 
before to set this saturation value to be of order unity, and therefore must reintroduce 
an explicit coefficient a in front of the Duffing term, which can no longer be scaled to 
unity. One can implement such an effect with an electronic feedback loop, sensing each 



oscillator velocity and driving the oscillator with an appropriate phase (Feng et al.. 



2008). The first three terms of eqn (11.5) comprise a so-called van der Pohl oscillator. 



Note that in anticipation of our study of synchronization of coupled oscillators in 



Sec. 11.4 below, we have assumed that the uncoupled oscillators can generally have 
non- identical linear frequencies aj„. 

The equations of motion for particular experimental implementations might have 
different terms, although we expect all will have positive or negative Duffing terms; 
linear and possibly also nonlinear damping; linear and possibly also nonlinear coupling, 
which may be either reactive or dissipative; and some source of energy to sustain the 
oscillations. In many cases, although not always, once we transform to the reduced 
description describing the slow modulation of the modes (see below), the differences 
between these different models will not lead to qualitatively new effects. 
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11.3 Discrete Amplitude Equations: 

Example I - Collective response to parametric excitation 

11.3.1 Deriving the equations 



As in the case of a single resonator in eqn (11.2), we suppose Q is large and take 
e ~ Q^^ as a small expansion parameter. Again, we take H = eh, and in addition also 
take D = ed so that the width of the frequency band of normal modes is also small. 
This is not quite how ^Lifshitz and Cross (2003| ) treated the coupling, but it is simpler 
yet equivalent up to the order of the expansion in e that we require. The equations of 



motion (11.4) then become 



+ eu„ + (1 + eh cos [(2 + eilp) t]) u„ + ^ed{u„ 



2u„ 



Un-l) 



■n + Wn'^n = 0. 

(11-6) 

We expand u„(i) as a sum of standing wave modes with slowly varying amplitudes. 
The nature of the standing wave modes will depend on the conditions at the ends of 
the array of resonators. In the experiment of Buks and Roukes (2002[ | there where 
N mobile beams with a number of identical immobilized beams at each end. These 
conditions can be implemented in a nearest neighbor model by taking two additional 
resonators, mq and u^+i and assuming 



Mo 



un+1 = 0. 



The standing wave modes are then 

Un = sin(ng„i) with g„ 



N +1' 



TO = 1 . . . iV. 



(11.7) 



fll. 



On the other hand, for an array of N resonators with free ends there is no force from 
outside the array. For the nearest neighbor model this can be imposed again by taking 
two additional resonators, but now with the conditions 



uq = ui; 



The standing wave modes are now 

ii„ = cos [(n - i) g„] with 



q-n 



N 



m = Q...N -I. 



For our illustration we will take eqns (11.7 
To treat the equations of motion 



11.8) 



(11.9) 



(11.10) 



11.6) analytically, we use secular perturbation 



theory combined with a multiple scales analysis, taking advantage of the natural sep- 
aration of time scales that occurs in our physical system — the fast oscillations of the 
resonators at half the drive frequency are characterized by the fast time variable 
whereas the slow variation of the amplitudes of these oscillations is associated with 
transient times, characterized by the damping rate , or e, giving rise to a well- 
separated slow time variable T = et. This approach is used throughout this review, 
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and was described in great detail in our previous review (Lifshitz and Cross, 2008 



Lifshitz and Cross, 2010). Thus, we introduce the ansatz 



N 



Mt) = ^^'^\ (An(T)sin(ng„)e** + c.c.)+e3/2uW(t) + ...^ n=l...N, (11.11) 

m— 1 

where c.c. stands for the complex conjugate. The lowest order contribution to this 
solution is based on the normal mode solutions (11.81 of the linear equations of motion, 
allowing the complex mode amplitudes Am{T) to vary slowly in time (as in a rotating 
frame in the complex plane), due to the effect of all the other terms in the equation. 
As we shall immediately see, the slow temporal variation of A^iT) also allows us to 
ensure that the pcrturbative correction Un\t), as well as all higher-order corrections 
to the solution ( 11.11 ), do not diverge as they do if one uses naive perturbation theory. 



Using the relation 



Ar. 



dAn 



dA„ 



(11.12) 



and denoting a time derivative with respect to the slow time T by a prime, we sub- 
stitute the trial solution (11.11) into the equations of motion (11.6) term by term. Up 
to order e^/^ we have. 



^ sin(ng„0 {[-A„, + 2ieA'Je^' + c.c.) + t^'^il^^Xt), 

m 



,e'* -I- c.c. 



(11.13a) 
(11.13b) 



e-[Un+l - 2Un + Un-l) = ~e ' - > 2 



sm 



2 f 1m 



) sin(ng„J (A,„e'* + c.c.) , (11.13c) 



e^/^^ X! sin(n(jfe) s'm{nqi) (A^e** + c.c.) (AfcC** + c.c.) {Aic'* + c.c.) 



3/2 



^ X! {sin[n(-gj + qu + qi)] + sm[n{qj - qk + qi)] + sm[n{qj + qu - qi)] 



32 



in[n{q, + qu + qi)]} {A.AkAic^'' + 3A,AkA*e'' + c.c.} , (11.13d) 



and 



rjUn^Un = e^^^T^I X! {sin[n(-(7j + qu + qi)] + sin[n(gj - qu + qi)] + su\[n{qj + qu - qi)] 



- s\n[n{qj + qu + qi)]} {Ajc'* + c.c.) (^fcc'' -|- c.c.) (iA,e" + c.c.) . (11.13e) 
The order e^/^ terms cancel, and at order e^^^ we get N equations of the form 

(m*'' secular term) e** -|- other terms, (11.14) 



where the left-hand sides are uncoupled linear harmonic resonators, with a frequency 
unity. On the right-hand sides we have N secular terms which a ct to drive the res- 
onators Un^ at their resonance frequencies. As Lifshitz and Cross (2008 2010) did for 
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all their single resonator examples, here too we require that all the secular terms vanish 
so that the Un"^ remain finite. This is the necessary solvability condition, required to 
obtain equations for the slowly varying amplitudes AmiT). To extract the equation for 
TO*'' amplitude Ajn{T) we make use of the orthogonality of the modes, multiplying 

th 



the 



all the terms by sin(7ig,„) and summing over n. We find that the coefficient of the to 
secular term, which is required to vanish, is given by 



-2i- 



It' 



iAm+idsin^ 



A„ 



1 



hA* 



JQpT 



IT] X ^ 



16 



A,Au.A1^ 



(1) _ 



(11-15) 

where we have used the A function introduced by Lifshitz and Cross (2003 ) , defined 
in terms of Kronecker deltas as 

.(1) 



^ j+k — l.m 



5j-k+l,-rn 
^j + k-l,-m 



^j-k+l,2iN+l)-m 
5j + fc-i,2(Ar+l)-m 



(11.16) 



— Sj+k+l,m + Sj+k+l,2{N+l)-m ^ Sj + k+l,2{N+l)+m, 



and have exploited the fact that it is invariant under any permutation of the indices j, 
fc, and I. The A function ensures the conservation of lattice momentum — the conser- 
vation of momentum to within the non-uniqueness of the specification of the normal 
modes due to the fact that sin(ng,„) — sin{nq2k(N+i)±m) for ^-i^y integer k. The first 
Kronecker delta in each line is a condition of direct momentum conservation, and 
the other two are the so-called umklapp conditions where only lattice momentum is 
conserved. 



As for the single resonator (Lifshitz and Cross, 2008), we again try a steady-state 
solution, this time of the form 



A^iT) 



(11.17) 



so that the solutions to the equations of motion (11.6), after substitution of (11.17) 
into ( 11.11[ ), become 



+ C.C. ] +0{e'^') 



(11.18) 



where all modes are oscillating at half the parametric excitation frequency, up 
2 + e^lp. 



Substituting the steady state solution (11.17) into the equations (11.15) for the 



time- varying amplitudes Am{T), we obtain the equations for the time-independent 
complex amplitudes a™ 



2c? sin 



. 2 f 1m \ 



h ^ S + ir]^^ (1) 
2 16~ ^ 

3:k.l 







(11.19) 



Note that the first two terms on the left-hand side indicate that the linear resonance 
frequency is not obtained for fip — 0, but rather for VLp + 2dsiT? {qm/2) = 0. In terms 
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of the unsealed parameters, this imphes that the resonance frequency of the m*'' mode 
is ujm = 1 — -Dsin^ (qm/^), which to within a correction of O(e^) is the same as the 
expected dispersion relation 



^1-2D sin' 



m 



(11.20) 



Equations (11.151 and (11.19) are the main result of the calculation. We have 



managed to replace N coupled differential equations (11.4) for the resonator coordi- 



nates Un{t) by N coupled differential equations (11.15) for the slowly varying mode 
amplitudes A,„(T), and then by N coupled algebraic equations (11.19) for the time- 
independent mode amplitudes a^. All that remains, in order to obtain the overall 
collective response of the array as a function of the parameters of the original equa- 



tions of motion (11.4), is to solve these coupled algebraic equations. 



11.3.2 Analyzing and solving the equations 

A number of simple results can immediately be stated. First, one can easily verify 
that for a single resonator (TV = j = k = l = m = l), the general equation (11.19) 
reduces to the single-resonator equation treated by Lifshitz and Cross ( 2008[ 2010), as 
Aiii;i = 4. Next, one can also see that the trivial solution, Um = for all to, always 



satisfies the equations, though, as Lifshitz and Cross (2008 20101 showed in the case 



of a single resonator, it is not always a stable solution. Finally, one can also verify that 
whenever for a given to, ^^rnm-j — ^ j fhexi a single-mode solution exists 

with am 7^ and aj = for all j ^ to. These single-mode solutions have the same type 

of elliptical shape of the single-resonator solution. Note that generically AmLmim = 3, 
except when umklapp conditions are satisfied. 




Fig. 11.1 Response intensity of two resonators as a function of frequency Q,p, for a particular 
choice of the equation parameters, (a) shows |ai|^, and (b) shows |a2p, with solid curves 
indicating stable solutions and dashed curves indicating unstable solutions, (c) Comparison 
of stable solutions, obtained algebraically (small circles), with a numerical integration of the 
equations of motion (11.41 (solid curve - frequency swept up; dashed curve - frequency swept 
down) showing hysteresis in the response. Plotted is the averaged response intensity, defined 
in the text. In all figures, the two elliptical single-mode solution branches are labeled Si and 
S2, and the two double-mode solution branches are labeled D\ and D2- From Lifshitz and 
Cross (2003). Copyright (2003) American Physical Society. 
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Additional solutions, involving more than a single mode, exist in general but are 
hard to obtain analytically. Lifshitz and Cross (2003[ ) calculated these multi-mode solu- 



tions explicitly for the case of two and three resonators, for the model they considered, 
by finding the roots of the coupled algebraic equations numerically. We present some of 
their results to illustrate the type of behavior that occurs, although the precise details 
will be slightly different in the model used here. In Fig. |11.1| we show the solutions 
for the response intensity of two resonators as a function of frequency, for a particular 
choice of the equation parameters. Figure ll.l[ a) shows the square of the amplitude 



of the symmetric mode oi, whereas Fig. ll.l[ b) shows the square of the amplitude of 



the antisymmetric mode 02. Solid curves indicate stable solutions and dashed curves 
indicate unstable solutions. Two elliptical single-mode solution branches, similar to 
the response of a single resonator are easily spotted. These branches are labeled by 



and ^2. Lifshitz and Cross (2003) give the analytical expressions for these two solution 
branches. In addition, there are two double-mode solution branches, labeled Di and 
D2, involving the excitation of both modes simultaneously. Note that the two branches 
of double-mode solutions intersect at a point where they switch their stability. 

With two resonators there are regions in frequency where three stable solutions can 
exist. If all the stable solution branches are accessible experimentally then the observed 
effects of hysteresis might be more complex than in the simple case of a single resonator. 
This is demonstrated in Fig. |ll.l| ^c) where the algebraic solutions are compared with a 



numerical integration of the differential equations of motion (11.4) for two resonators. 
The response intensity, plotted here, is the time and space averages of the square of 
the resonator displacements {{uf) + {u'^))/2, where the angular brackets denote time 
average. A solid curve shows the response intensity for an upward quasistatic frequency 
sweep, and a dashed curve shows the response intensity for a downward sweep. Small 
circles show the response intensity, as calculated for the stable regions of the four 
algebraic solution branches shown in Figs. [TLl^ a) and (b), demonstrating the great 
utility of the slow amplitude equations. With the analytical solution in the background, 
one can easily understand all the discontinuous jumps, as well as the hysteresis effects, 
that are obtained in the numerical solution of the equations of motion. Note that the 
5*1 branch is missed in the upward frequency sweep and is only accessed by the system 
in the downward sweep. One could trace the whole stable region of the Si branch by 
changing the sweep direction after jumping onto the branch, thereby climbing all the 
way up to the end of the Si branch. These kinds of changes in the direction of the 
quasistatic sweep whenever one jumps onto a new branch are essential if one wants 
to trace out as much of the solution as possible — whether in real experiments or in 
numerical simulations. 

11.3.3 Brief survey of applications 



Discrete amplitude equations like the ones derived here ( 11.19 ) are useful mainly when 
studying small arrays. Nevertheless, the insight gained by studying small arrays, of 
even two or three resonators, provides better understanding of the dynamics of large 
arrays, which can be studied directly by numerically integrating the starting equations 



of motion ( 11.4). Indeed, all the features of the original experiment of Buks and Roukes 



72002) were qualitatively reproduced by ILifshitz and Cross (20031) by numerically 
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integrating their original equations of motion (11.3). But, it was their analytical study 



of small arrays, that allowed them to provide an explanation for the observed features: 
(1) The response of the system at frequencies above the top edge of the band was 
attributed to the positive frequency pulling coming from the Duffing nonlinearity; (2) 
The fact that only a few features are observed in a monotonic quasistatic frequency 
scan, rather than A'' resonance peaks for the N normal modes, was explained by the fact 
that a solution branch is followed quasistatically as long as it is stable, often skipping 
many other solutions that are simultaneously stable, as demonstrated above with 2 
resonators; and (3) The abrupt jumps in the response were identified as stemming 
from bifurcation points where a certain solution branch ends, as in the saddle-node 



bifurcation at the end of the Di branch in Fig. 11.1 or simply loses its stability, as for 
the two S branches in Fig. |11.1[ in either case requiring the system to switch abruptly 
to a different branch. 

The curious reader is encouraged to consult additional articles, where the methods 
presented in this section were used to study advanced features in the dynamics of 



small numbers of coupled resonators. Karabalin et al. (2009 1 used discrete amplitude 
equations for the two normal modes of a pair of resonators, similar to eqn (11.151, to 



assist in their numerical modeling of period doubling and a transition to chaos, which 



they observed experimentally. Kenig et al. (2011 1 used discrete amplitude equations to 
identify homoclinic orbits in the slow dynamics and assess the possibility of obtaining 



chaotic dynamics via the Melnikov approach. Finally, Karabalin et al. (2011) demon- 
strated the use of a pair of parametrically driven resonators as a novel amplifier, whose 
operation is based on very sensitive control of the bifurcation diagram of the response 
of two resonators, via an input signal that is fed into the coupling D between the 



resonators. The Supplementary Material of Karabalin et al. (2011 ) provides a detailed 



analysis of the operation of this so-called Bifurcation- Topology Amplifier, using a set 
of discrete amplitude equations like the ones developed here. 



11.4 Discrete Amplitude Equations: 

Example II - Synchronization of nonlinear oscillators 

11.4.1 Deriving the equations 

Although synchronization is often put forward as an example of the importance of 
understanding nonlinear phenomena, the intuition for it, and indeed the subsequent 
mathematical discussion, often reduces to simple linear ideas. For example, the famous 



example of Huygens's clocks (Bennett et al., 2002) can be understood in terms of a 



linear coupling of the two pendulums through the common mounting support. It is 
then the larger damping of the symmetric mode (coming from the larger, dissipative 
motion of the common support) compared with the antisymmetric mode that leads, 
at long times, to a synchronized state of the two pendulums oscillating in antiphase. 
The nonlinearity in the system is simply present in the individual motion of each pen- 
dulum; specifically in the mechanism to sustain the oscillations. Without the drive, 
the oscillators would still become synchronized through the faster decay of the even 
mode, albeit in a slowly decaying state. Rather than this mode-dependent dissipation 
mechanism, one might expect synchronization to arise from the intrinsically nonlinear 
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effect of the frequency pulling of one oscillator by another. Furthermore, the model 
describing the two Huygens pendulums, as well as most other models used to show 
synchronization, has dissipative coupling between the oscillators. In contrast, many 
physical situations have mainly reactive coupling. Consequently, Cross et al. (2004 
[2OO6 ) proposed and analyzed a model for synchronization, given by eqn (11.5), in- 



volving reactive coupling between the oscillators, which then leads to synchronization 
through nonlinear frequency pulling. 



We follow Cross et al. ( 2004 2006 ) and consider the system of oscillators defined 



by eqn (11.51, assuming that the linear frequencies of the oscillators are distributed 



near unity such that 



= 1 + A„, 



with 



|A„| « 1. 



(11.21) 



This allows us to study the situation in which the equations of motion are dominated 
by the terms describing simple harmonic oscillators at frequency one, and the time 
dependence remains close to e^**. The interesting dynamics should then be captured 
by a discrete set of coupled amplitude equations for the deviations of the individual 
oscillators from simple harmonic oscillation at frequency 1. To that end we assume that 



all corrections in eqn (11.51 to a set of uncoupled harmonic oscillators of frequency 1 



are small. To formalize this smallness we again use the damping term to define a small 
parameter e = j^, and take A„ = e^„, a = ea/3, D = e/3. The oscillating displacement 
is then written as a slow modulation of oscillations at frequency one, plus corrections 



u„{t) = [A,,{T)e^' + c.c] + euil\t) + 



(11.22) 



with T = a slow time scale as before. As always, the slow variation of An{T) gives 
us the extra freedom to eliminate secular terms and ensure that the perturbative 
correction Un\t), as well as all higher-order corrections to the linear response, do 
not diverge. Note that our decision to scale eqn (11.5) by setting the van der Pohl 
term such that the nonlinear saturation of the oscillations occurs at u„ = 0(1), has 



affected the scaling of our trial solution ( 11.22 ), whose leading term is indeed of order 



1. Compare this with the trial solution of the previous section, given by eqn (11.11), 
whose leading term is of 0(\/e). 



Using the relation (11.12), again denoting a time derivative with respect to the 



slow time T by a prime, we calculate the time derivatives of the trial solution ( 11.22 ) 

iin = {[iA. 



eA' 



y*+c.c.) 



eu, 

it 



(1) 



([-A„ + 2ieA'^ + €^A';y + c.c.) + 



eu. 



(1) 



it) 



(11.23a) 
(11.23b) 



Substituting these expressions back into the scaled equation of motion 



■tin + (1 + e6n) U„ 



[(1 



\aul 



2Ur. 



,_i)] =0, (11.24) 



and picking out all terms of order e, we get the following equation for the first pertur- 
bative correction 



,(1) 



5nAn - (2iA'„e" + c.c) + (iA„e** + c.c) 
\a (A„e** + c.c.f + \p [{An+i - 2A„ + 



1 



i)e'* 



+ c.c. 



+ c.c. 



(11.25) 
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Terms varying as e on the right-hand side of Eq. (11.251 contribute a finite 



response to Un \ but the collection of terms proportional to e'* — the secular terms — act 
like a force driving the simple harmonic oscillator on the left-hand side at its resonance 
frequency. The sum of al l these secular terms must vanish so that the perturbative 
correction Un\t) in eqn (11.22) will not diverge. This provides us with a solvability 
condition that leads to an equation for determining the slowly varying amplitudes 

An{T) 

2^ = (1 + iSn)An - (1 - za) |A„|' An + - 2An + . (11.26) 

With a rescaling of time t = T/2 = et/2 and a slight rearrangement of terms, 
eqn (11.261 reduces to the form obtained by Cross et al. (2004 20061, 



^ = i{5n + a |A„|')A„ + (1 - |A„|')A„ + {An+i - 2A„ + A„_i) 



(11.27) 



The first term on the right-hand side shows the ability of the n*'' oscillator to 
shift its frequency by an amount a |A„|^; the second term shows the tendency of the 
oscillators to increase their amplitude as long as |A„|^ < 1; and the third term is 
the reactive coupling between nearest neighbors. If this nearest-neighbor coupling is 
generalized to allow also dispersive interaction and replaced by an all-to-all or mean- 
field coupling, convenient for theoretical analysis, we obtain the final form of the model 



studied by Cross et al. (2004 2006) 




(11.28) 



except for the fact that we use the opposite sign convention for the Duffing parameter. 
Thus, in our current discussion a positive (negative) value of a implies a stiffening 
(softening) Duffing nonlinearity. The relative natural frequency i5„ of each oscillator is 
chosen from a specified distribution g{5)^ whose width is denoted by w. 

When only nonlinear saturation and dissipative coupling are present [a = (3 = 
0, if ^ 0) eqn (flT^ reduces to 



— — = (Z(5„ + 1 

dr 



\An\ )j4„ 



m— 1 



(11.29) 



which has been analyzed by Matthews et al. (1991) for general w and K. 
11.4.2 Analyzing and solving the equations 

The complex number An, representing the amplitude r„ and phase On of the n*'' 
oscillator, An = r„e'^", suggests the introduction of a complex order parameter 4' to 
measure the coherence of the oscillations 
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* = i? e'® = 



1 ^ 
N ^ 



(11.30) 



A nonzero value of the order parameter i? > may be taken as the definition of a 
synchronized state. 

The general amplitude-phase model reduces to familiar phase only models of syn- 
chronization in certain limits. If the width w of the distribution g{S) is narrow, so 
that the time evolution of the magnitudes r„ = |A„| is fast compared with that of the 
phase dispersion, and the coupling constants K, (} are small, r„ rapidly relaxes to a 
value close to unity 



K 

N 



N 

El 



/3 



N 



m—1 



sm 



(dm — Sn), 



(11.31) 



and the only remaining dynamical variable for each oscillator is its phase 9„. Equa- 



tion (11.281 can then be reduced to 



9„ = Sn + a + 



N 



N 

E 

m—1 



sm 



N 



N 

E 

m—1 



COS 



0„)~1]. (11.32) 



For the case of purely dissipative coupling a = /3 = 0, if 7^ 0, or reactive coupling with 
strong frequency pulling K — Q,a, (3 ^ Q, \a\ ^ 1, the last term on the right hand side 
of eqn (11.32) can be neglected, and the equation reduces to a simple form (Winfree, 
1967[ Kuramoto, 1975), known as the Kuramoto model (ignoring the unimportant 
constant term a and writing the effective coupling constant in either case simply as a 
K) 

N 

E 

m—1 



= 5n 



K 

N 



sm 



(11.33) 



that has been the subject of numerous studies (Acebron et al., 2005). In the absence of 



coupling each oscillator in this model would simply advance at a rate that is constant 
in time, but with some dispersion of frequencies over the different elements. 

Identifying the imaginary part of 
while recalling that r„ = 
mean-field expression 

9^^6n + KRam{e-e„). (11.34) 



in the sum appearing in eqn (11.33) 
for the Kuramoto model 



-yields a particularly simple 



Thus the behavior of each oscillator is given by its tendency to lock to the phase of the 
order parameter. The term KRsm{ld — On) acts as a locking force, and locking occurs 
for all oscillators with frequencies satisfying |(5„| < KR, with the locked oscillator 
phase given by -|- sm~^{5n/ KR). The magnitude R of the order parameter must 



then be determined self-consistently via eqn (11.30) 



Equation (11.33) is known to show rich behavior, including, in the large N limit, 
a sharp synchronization transition at some value of the coupling constant K = Kc 



(Kuramoto, 1975), which depends on the frequency distribution g{5) of the uncoupled 



14 Collective Dynamics in Arrays of Coupled Nonlinear Resonators 



oscillators. The transition is from an unsynchronized state with = in which the 
oscillators run at their individual frequencies, to a synchronized state with 7^ in 
which a finite fraction of the oscillators lock to a single frequency. The transition at 
Kc has many of the features of a second order phase transition, with universal power 
laws and critical slowing down (Kuramoto, 1975), as well as a diverging response to 



an applied force (Sakaguchi, 19881 



The last term in eqn (11.321 may lead to important qualitative effects even if the 



coefficient is not very large. For example, in the absence of this term the coupling terms 
cancel when summing over all the oscillators in the system, so that the frequency of a 
synchronized state is simply related to the mean frequency of the oscillators. This is no 
longer the case for the general equation. In the case of short range, rather than all to 
all coupling, the cos{9m — dn) term profoundly changes the nature of the synchronized 



state to one of propagating waves (Sakaguchi et al., 1988 Blasius and Tonjes, 2005) 



11.4.3 Brief survey of applications 

The synchronization of oscillators with reactive all-to-all coupling and nonlinear fre- 
quency pulling, described by of eqn (11.28) with K = 0, was analyzed by Cross et 
al. (2004 ,2006) for several different frequency distributions g{5) (Lorentzian, top-hat. 



and triangular). Here we briefiy review the results for a triangular distribution with 
width w = 2, and refer the reader to the original work for more details. Such a width 
is not small compared with the relaxation rate of the magnitude variables, and so 
the behavior is richer than in the weak randomness limit described by the Kuramoto 
model. The stability diagram of the variety of states found as a and (3 are varied is 
shown in Fig. [TL2| ;d). The resuhs are shown for a/3 < 0, noting that for a symmet- 
ric distribution of frequencies the results are the same if both signs of a and /3 are 



changed. These same results were presented by Cross et al. { 2004 2006 ) for their case 



aP > 0, as they were using the opposite sign convention for a. 

Certain results can be obtained analytically, in particular the instability from the 
unsynchronized state to a synchronized state with nonzero order parameter i? > 0, 
and the instability from the fully locked state (all oscillators locked to evolve at the 
same frequency) to a synchronized state with only partial locking. Other results are 
obtained numerically by performing sweeps of /? at fixed values of a. Results for the 
time averaged magnitude of the order parameter {R)t for three values of a are shown 
in Fig. 11.2[ a)-(c). Both upward and downward quasistatic sweeps of the reactive 
coupling strength /3 are used to uncover hysteresis in the transitions. 

Many novel features are apparent in these results. For example, the unsynchronized 
state is stable for both small and large values of the coupling strength /3, so that for 
fixed a there are two values of /3 at which the unsynchronized state passes from stable 
to unstable. At large values of /3 a large order-parameter synchronized state is also 
stable, which becomes the fully locked state at large enough /3. The transition from the 
unsynchronized state to the synchronized state may be continuous, passing through a 
supercritical bifurcation, or discontinuous, passing through a subcritical bifurcation. 
Hysteresis is apparent in the latter case, owing to the bistability of both the synchro- 



nized and the unsynchronized states. This is demonstrated explicitly in Fig. 11.2 ^c). 
More surprisingly, one also observes the multistability of different synchronized states. 
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Fig. 11.2 (a)-(c) Simulations of 1000 oscillators having a triangular frequency distribution 
with width w — 2. The time-averaged order parameter magnitude {R)t is plotted for both 
upward and downward sweeps of /3 at fixed a: (a) a = 0.0; (b) a = —0.4; and (c) a = —0.9. 
The same results would be obtained by switching the signs of both a and p. (d) Stability 
diagram for the same triangular distribution, in the quadrant of the a — /3 plane with a/3 < 0. 
Solid and dashed lines show analytical results of the linear stability of the unsynchronized 
state. Numerics show the bifurcations are supercritical along the solid portions and subcritical 
along the dashed portion. Dotted line is the linear stability boundary of the fully locked state. 
Dash-dotted lines are saddle-node bifurcations observed in numerical simulations. States are: 
U - unsynchronized; 5*1 , 52 synchronized with small and large amplitude respectively; L fully 
locked. From Cross et al. (2006). Copyright (2006) American Physical Society. 

Over some parameter ranges a small order-parameter synchronized state may coex- 
ist with the large order-parameter synchronized state, as observed, for example, in 
Fig. |11.2| c) between about j3 = 3.0 and (3 = 3.7. This small order-parameter synchro- 
nized state has the novel property that the order parameter is nonzero R > Q, but 
there is no oscillator locked in frequency to the frequency of the order parameter or 
of other oscillators — this is a synchronized state i? > with no frequency locking. 

The rich synchronization behavior displayed by this stability diagram opens up 
many possibilities for applications, as well as suggesting difhculties that must be over- 
come, for example when there exists a multistability of different dynamical states. 



11.5 Continuous Amplitude Equations: 

Example III - Nonlinear competition between extended modes 

11.5.1 Derivation of the BCL amplitude equation 

We wish to investigate the sequence of single mode standing wave patterns to be ex- 
pected in parametrically driven resonator arrays, in cases where many such modes are 
simultaneously stable, when the strength of the drive is varied. Although the quantita- 
tive analysis could be done directly from the basic equations of motion for the coupled 
resonators, it is advantageous to formulate the analysis in terms of a continuous am- 
plitude equation — that which was developed by Bromberg, Cross and Lifshitz (2006), 
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henceforth referred to as the BCL equation. This aUows us to display the range of 
stable patterns on a reduced stability diagram involving just two dimensionless vari- 
ables (a scaled measure of the driving strength, and a scaled mode wave number), 
so that it is easy to deduce the general qualitative behavior upon variation of the 
parameters. The specific quantitative behavior for a physical system is also easy to 
obtain by evaluating the corresponding scaled quantities. A change of pattern occurs 
when parameters vary so that the mode moves outside of the region of stable patterns 
on this diagram, and the new pattern is predicted by analyzing the result of the in- 
stability using the BCL equation. This type of approach was used in other pattern 



forming systems (Kramer et al., 1988). A novel feature of the present system is that 



the difference in the instabilities encountered on increasing and decreasing the (scaled) 
driving strength leads to the prediction of quite different-sized mode jumps for the up 
and down sweeps. 

We follow the treatment of BCL in deriving their amplitude equation, but instead 



of starting with the original equations of motion (11.3) derived by Lifshitz and Cross 
[(2003^ , we start with the simpler equations of motion ( 11.4 ). This leads to a somewhat 
simplified derivation, which eventually yields the same amplitude equation to describe 
the slow dynamics of the system of resonators. We perform the same scaling of the 
equation parameters as we did in Sec. 11.3 with one difference — in anticipation of 



treating extremely large arrays, with thousands or more normal modes of vibration, 
we do not wish to assume that the width of the frequency band is small. We therefore 
do not replace D with ed as before. This will also allow us to obtain the exact dispersion 
relation (11.20) at the linear step, i.e. at order y^. Our starting point is therefore the 



set of coupled equations 

Un + eUn + (1 — ehcOs2ujpt) -U„ + ^D{Un+l — 2Un + U,i-l) + ufi + "qu^Un = 0, (11.35) 

where we have taken a negative sign for the parametric driving term to be consistent 
with the sign used by BCL, thus merely shifting the phase of the drive by tt relative 



to eqn (11.41 



Amplitude Equations for Counter Propagating Waves. In order to treat this system 
of equations analytically, beyond the treatment described earlier in Sec. 11. 3[ we in- 
troduce a continuous displacement field u{x,t), keeping in mind that only for integral 
values X = n oi the spatial coordinate does it actually correspond to the displacements 
u{n, t) = Un{t) of the discrete set of resonators in the array. We introduce slow spatial 
and temporal scales, X = ex and T — et, upon which the dynamics of the envelope 
function occurs, and expand the displacement field in terms of e. 



m(x, t) 



= .1/2 



[(A+(A,T)e^'«''^ Al(X,r)e*'?''^) e*"''* + c.c." 



+ e'/^u('Hx,t,X,T) + 



(11.36) 



where the asterisk and c.c. stand for the complex conjugate, and qp and ujp are related 
through the dispersion relation. 



(11.37) 
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Note that the response to lowest order in e is expressed in terms of two counter- 
propagating waves with complex amplitudes and A_ , which is a typical ansatz for 
parametrically excited continuous systems ( [Cross and Hohenberg, 1993 ) . We substitute 
the ansatz (11.361 into the equations of motion ( [11.35 ) term by term. Again, using 
eqn ( [11.12[ ) in addition to expanding A±{X + e, T) ~ A±{X, T) + edA±{X, T)/dX we 
obtain up to order e'^/^, 



.1/2 



+ c.c. + e 



3/2 



dA+ 



+ ~ujt,A*_ + 2iujpe- 



dA*_ 
W 



(11.38a) 



Un±l 



.1/2 



A. ±e 



dA, 



dX 

+ e^/^u'^'^\x±l,t,X, T), 



dA*_ 
~dX 



Jqp{x±l) 



e*"p* + c.c. 



(11.38b) 



-D{u. 



n+l 



2u„ + u„_i) = -e^/^2Dsm^{qp/2) (A+e-'"^"^ + A*_e"i^^) e*""* 



e^/^iDsin{qp) 
,D 



fdA+^_,^^.^ dA 



iqpX \ iujpt 



+ C.C. 



,3/2: 



\dX " dX 
{x+l,t,X,T) - 2u(i) (a;, t, X, T) + u^^^ (a: - 1, <, X, T) 



e/icos(2wpi)u„ 6^/2- (A^e^"'"^ + ^;e"'''^) e'""* + 0{c'^'^^') + c.c, (11.38d) 



■3 - + 2|A_|2) A+c-^i^"" + (2|A+|2 + Ale'*"^] 



(11.38c) 
(11.38d) 
(11.38e) 



+ c.c, 



(11.38f) 



and 



+ O (e^^-'p*, e^^'"^) + c.c, (11.38g) 

where 0{e^^"p* , e^^'^p^) are fast oscillating terms proportional to e*^"^"* or e*'^*"'^ that 
do not enter the dynamics at the lowest order in e beca use they are nonsecular. 

At the order of ei/^, the eq uations of motion ( 11.35 ) are satisfied trivially, yielding 
the dispersion relation (11.37) mentioned earlier. At the order of e"^/^ on the other 
hand, we again obtain secular terms, and must apply a solvability condition, which 
requires that all terms proportional to g^i'^pt^qpx) jj^^gf^ vanish. As a result, we obtain 
the two coupled amplitude equations. 



dA± dA± 

± Vo 



dT 



dX 



1 ih 1 

"2^±^4^^^" 2 



3i 



ryT— (|A±|2 + 2|A^|2)^±, (11.39) 



where the upper signs (lower signs) give the equation for ( A_ ) , from the restriction 
on the terms proportional to e"^p*~*9p^ (^^u^pt+iqpx-^^ g^j^j where 
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duJ 
dq 



D sin(gp) 



(11.40) 



is the group velocity. A detailed derivation of the amplitude equations (11.391 can be 



found in the Masters thesis of Bromberg (2004). Similar equations were previously 
derived for describing Faraday waves (Ezerskii et al., 1986 Milner, 1991). 



Reduction to a Single Amplitude Equation. By linearizing eqns ( 11.39 ) about the zero 



solution (A^ — A- — 0) we find that the linear combination of the two amplitudes that 
first becomes unstable at he — 2LUp is i? oc (yl-|_ — iA-) — representing the emergence 
of a standing wave with a temporal phase of 7r/4 relative to the drive — while the 
orthogonal linear combination of the amplitudes decays exponentially and does not 
participate in the dynamics at onset. Thus, just above threshold we can reduce the 
description of the dynamics to a single amplitude B, where at a finite distance above 
threshold a band of unstable modes around qp can contribute to the spatial form of B. 
This is similar to the procedure introduced by Riecke (1990) for describing the onset 
of Faraday waves. 

To proceed with our multiple scales analysis, and obtain an equation describing the 
relevant slow dynamics of the new amplitude B, we need to identify another physically 
small parameter with which we can associate even slower spatial and temporal scales. 
We therefore assume that the coefficient of nonlinear damping rj is small, and define 
a second small parameter 5 = 77^ <C 1. We then define a reduced driving amplitude g 
with respect to the threshold he by letting (/i — he) /he = gS. A sequence of judicious 
arguments ( Bromberg, 2004} Bromberg et al, 2006[ ) then encourages us to scale the 



original amplitudes A± as S'^^^, making the ansatz that 



^1/4 



,(1) 



(^,r,e,f) 



«(2)(x,T,e>) 



where ^ = S^^'^X and f 



(11.41) 

ST are the new spatial and temporal scales respectively. 



We substitute the ansatz (11.41) into the coupled amplitude equations (11.39) 



and collect terms of different orders in 6. Again, to the lowest order of e xpansi on the 
equations are satisfied trivially. Collecting all terms of order S^^'^ in eqns ( 11.39 ) yields 



^(1) 

yd) 



-v„^+z^\B\'-B 



1 



where exactly at onset D is a linear operator given by the matrix 

'5t + Vgdx + \ 



Ot - Vgdx + \ 



(11.42) 



(11.43) 



The vector ( J^j); on the right hand side of eqn (11.42) is an eigenvector of D, with 
an eigenvalue —1. The solution of eqn (11.42) is therefore immediately given by 



,(1) 



2uj„ 



(11.44) 
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We substitute eqn (11.441 back into eqns (11.39), collect all the terms of order S^^^ 
and obtain 



(2) 
.(2) 



W 



dB 
~df 



,d^B 



vz^ + -B--\B\^B 

9 i^c2 2 2 



4IBI 



,dB 



B 



.dB* 



2u;r 



(11.45) 



The vector ^ ^ , on the right hand side of eqn ( 11.45 ), is an eigenvector of D with zero 



eigenvalue. Clearly, the left hand side of the equation cannot contain any component 
along the direction of such a zero eigenvector. Therefore, the expression within the 
square brackets is a secular term that must vanish. This provides us with the required 



solvability condition to proceed ( Cross and Hohenberg, 1993 1 . After applying one last 
set of rescaling transformations, 
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(11.46) 



we end up with the BCL amplitude equation, which is governed by a single parameter, 



dB 



d^B .2 



4\B\ 



,dB_ 



B 



,dB* 



(11.47) 



11.5.2 Analyzing and solving the equation 



The simplest nontrivial solutions of the BCL amplitude equation (11.47) are steady- 
state single-mode extended patterns, given by 



(11.48) 



with bk and cp both real, and where the boundary conditions u(0,t) = u{N + l,t) = 
constrain the phase ip to be tt/4 or 57r/4. In steady state, the relation between 



the magnitude bk and the wave number k is found by substituting eqn (11.48) into 



eqn (11.47), and setting the time derivative to zero to give 



bl = {k-l) + V(fc-l)2 + (.9-fc2) > 0, 



(11.49) 



along with a negative square-root branch which is always unstable against small per- 
turbations, as can be verified by the analysis below. 

Substituting the single- mode solution of eqn (11.48), with ip = tt/A, back into 



eqn (11.44) and eqn (11.41), and then into eqn (11.36), yields extended single-mode 
standing-wave parametric oscillations at half the drive frequency, whose explicit form 
is given by 



/ .-s i/2ri/4 4wpv/l + tan2(a) m ko\ 

i{x, t) ~ e ' ' — — Ofc sin[qmX) cos(7r/4 — ujpt — a), (11.50) 

3v 3 
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where we have defined 



tan(a)^<5i/2'^(52_fc) 



To satisfy the boundary conditions m(0,<) 
must be of the form 



+ 1 



(11.51) 

u{N + 1,<) = 0, the wave numbers q„i 

(11.52) 



where 



AQ 



N 



3Dsm{qp) TT 



£,51/2 



1 



(11.53) 



BCL showed that the first single-mode pattern to emerge as the zero-state becomes 
unstable is the one whose wave number g,„ is closest to the wave number qp that is 
determined by the drive frequency ojp through the dispersion relation (11.371. This 
determines the value of the scaled wave number in the single- mode solution ( 11.48 ) to 
be 

ko=(m~ 9p^^) AQjv, (11.54) 

where m is the integer closest to qp{N + l)/7r. Note that AQn tends to zero as the 
size N of the array of resonators tends to infinity. 



Linearization of the BCL amplitude equation ( 11.47) shows that the zero state with 
B{£,,t) = — which is a solution of eqn ( 11.47 ) for any value of g — is stable against the 
formation of single-mode patterns with wave number k as long as g < fc^ . The neutral 
stability curve g = k"^ is plotted as a dashed parabola in Fig. 11.3 Furthermore, 
for fc < 1 the bifurcation from the zero state to that of single-mode oscillations is 
supercritical, occurring on the neutral stability curve, while for A; > 1 it is subcritical, 
with a locus of saddle-node bifurcations located along the line g = 2k — 1 (shown in 
Fig. 11.3 as a solid green line), where the square root in eqn (11.49) is exactly zero. 

The stability of a single-mode solution (11.481 of wave number k against an Eck- 



haus transition to a different single-mode solution of wave number ± Q is found by 
performing a linear stability analysis of solutions of the form 

B(e,r) = bkc-'''^ + (/3+(r)e-^('=+'3)« +/3* (T)e-'('=-'3)«) , (II.55) 

with 1/3-1- 1 <^ 1. When the larger of the two eigenvalues describing the growth of such a 
perturbation becomes positive the single-mode solution of wave number k undergoes 
an Eckhaus instability with respect to single-mode solutions of wave numbers k ± Q. 

For an infinite number of oscillators the Eckhaus instability forms the upper bound- 
ary of the stability balloon of the single-mode solutions, and also the lower boundary 
for k < 5/2. For k > 5/2 the lower boundary is the saddle node bifurcation line. 
For a finite number of oscillators, restricting Q to be an integer multiple of AQat in 
eqn (11.55) slightly shifts the Eckhaus instability lines. The upper Eckhaus bound- 



ary is shifted to larger values of g. The nature of the lower instability boundary now 
depends on the number of resonators in the array through AQn, as well as on the 
wave number fc. For fc < 1 the lower boundary will be the Eckhaus instability curve 
if |fc| > and the neutral stability curve otherwise. Because the only wave 
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number to satisfy \k\ < AQn/2 is ko, given by eqn (11.54), upon decreasing g the 
fcg solution undergoes a continuous transition to the zero state. For A: > 1 the lower 
boundary will be the Eckhaus instability curve if 1 < fc < (5 — 3{AQn/2)'^)/2, and 
the line of saddle node bifurcations otherwise. For AQ^ > 2 there is no portion of 
Eckhaus instability on the lower boundary, which is the neutral stability curve if fc < 1 
and the saddle node bifurcation curve if fc > 1. These stability boundaries are shown 
in Fig. |11.3| for an infinite system and for a system of = 92 resonators. Further 
details can be found in Bromberg et al. (2006| ) and Kenig et al. (2009a). 



11.5.3 Brief survey of applications 



Kenig et al. {2009a) used the BCL amplitude equation (11.47) to study a number of 
collective dynamical effects in one-dimensional arrays of coupled nonlinear resonators. 
The common thread linking these effects is the so-called question of pattern selection 




Fig. 11.3 (Color) Stability boundaries of the single-mode solution (11.481 of the BCL 
amplitude equation (11.471 in the g vs. k plane. Dashed line: neutral stability curve g — k^. 



Dotted line: stability boundary of the single-mode solution (11.48 1 for a continuous spectrum 
(AQjv — >■ 0). Solid lines: stability boundary of the single-mode solution for N — 92 and the 
parameters D = 0.25, Qp = 737r/101, and e = 5 = 0.01 (giving fco ==; -0.81 and AQjv ~ 3.70). 



Black line: the value of g for which perturbations of the form given by eqn (11.551 start to 
grow. Red line: the lower bound for k < 1, g = k^ . Green line: the lower bound for fc > 1, 
the locus of saddle-node bifurcations g — 2k — 1. Vertical and horizontal arrows mark the 
secondary instability transitions that are expected upon quasistatic sweeps of g. The blue 
upward-pointing arrows are for upward sweeps that undergo an Eckhaus instability, and 
the red downward-pointing arrows are for downward sweeps, of which the two with fc > 1 
experience a saddle-node bifurcation, and the one with fc < 1 goes through a continuous 
(supercritical) bifurcation. From Kenig et al. (2009a). Copyright (2009) American Physical 
Society. 
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( Cross and Greenside, 2009 ) — the nonlinear competition between different single-mode 



standing- wave patterns of the form of eqn (11.50), when many such solutions are 



simultaneously stable. This question becomes particularly interesting when the control 
parameter — in our case the drive amplitude g — is varied as a function of time, either 
quasistatically, abruptly, or in an intermediate ramp rate. In all such cases one is 
interested in predicting which of all stable patterns will be selected, as well as in 
the detailed understanding of the nature of the switching between patterns as their 
stability changes. 



The BCL amplitude equation allowed Kenig et al. (2009a) to map out the expected 
behavior of the resonators using universal stability diagrams, like the one shown in 



Fig. 11.3 Such a diagram immediately shows the type of instability that will be en- 
countered upon variation of the control parameter, and gives qualitative insights on 
the mode jumps to be expected. For example, for quasistatic parameter variations the 
jump in the mode number is always unity if the control parameter is increased so 
that the Eckhaus instability operates, but larger jumps are often seen if the control 
parameter is decreased so that a saddle-node bifurcation occurs. This is indicated by 
the dashed sideways arrows in Fig. |11.3 



It is instructive to describe how |Kenig et al. (2009a ) examined the process by which 
these two types of pattern switchings occur. To do so, we expand the general solution 
of the BCL amplitude equation in the linear modes of the array 



i3(e,r)=5]6„(T)< 



(11.56) 



where fc„ = fco + nAQjy, and fco is defined in eqn (11.54). Substituting a truncated 



mode expansion (11.56) containing a finite number of modes around feg into the BCL 



amplitude equation (11.47), allows us to replace this partial differential equation with 



a finite number of ordinary differential equations for the coupled mode amplitudes. 



dbn 
Or 



m.p ^ 

— ^ bmbibpbrb^_i_^_pj^_^,_^ 



AQn bmbph: 



m-\-p—n 



(11.57) 



To satisfy the boundary conditions, as mentioned above for the single-mode solu- 



tion (11.48), we take each mode amplitude to be zero at the boundaries by setting all 



the phases (p„ in Eq. (11.56) to 7r/4, and take the amplitudes &„ to be real, keeping 
in mind that they can be either positive or negative. Note that if all mode amplitudes 
except 6o are set to zero we obtain a single equation with n = m= p = l — r — 0, 



whose steady-state solution is the same as the single- mode solution of BCL (11.49) 



We take a closer look at the transient behavior during the first Eckhaus transition 
from the initial fco pattern to the fci pattern by plotting the time evolution of the four 



largest modes, as shown in Fig. 11.4 a). One can observe the decay of the unstable 



mode amplitude bo followed by the growth of bi to its steady-state value. One can also 
see that during the transient the amplitude of the unstable mode 5_i becomes non- 
zero. Its participation in the Eckhaus transition from the fco pattern to the fci pattern 
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is essential, as can be verified by considering these two modes alone in a truncated 
expansion. Limiting the expansion to bo and 61 suppresses the Eckhaus transition, 
and the ko pattern remains stable as g exceeds its expected value for the Eckhaus 
instability. The Eckhaus transition is observed only when the k-i mode is included as 
well, corresponding to the stability calculation, performed earlier for the state given 
by Eq. (ITr55|. 



One might naively expect that the same mechanism causes the transition from 
the ks pattern to the ki pattern at g = 19 through a double phase slip (with Q = 
2AQn), however, this is not the case. Fig. 11.4| [b) reveals the transient processes on 
a downward sweep of g just below the saddle node at g = 19. As g crosses the saddle- 
node bifurcation point, the amplitude 63 drops abruptly to zero. As can be seen from 



Eq. (11.571, in the zero-displacement state the linear growth rates of the solutions 
(11.56) are g — fc^, so the ko pattern has the largest possible growth rate and it out- 



grows the other modes until its amplitude approaches the steady state value (11.49) 



However, at this value of g the ko pattern is Eckhaus unstable with respect to the fci 
pattern — notice the characteristic evolution of the modes around r = 3 in Fig. 11.4[ b) 
corresponding to the Eckhaus instability [cf. around t = 25 in Fig. 11.4[ a)]. Thus the 
fci mode is ultimately the selected pattern. 

For more rapid increases in the control parameter larger jumps in the mode number 
may occur, and these were shown to be predicted simply from a linear stability analysis 
following the Eckhaus instability. Kenig et al. (2009a| ) showed that following an abrupt 



10 15 20 25 30 35 40 



(a) Eckhaus instability 



(b) Saddle-node bifurcation 



Fig. 11.4 Time evolution of the amplitudes of the four largest modes that participate in 
(a) the transition from the initial ko pattern to the fci pattern, when the value of the control 
parameter is changed from <; = 10 to (; = 11, causing the initial fco pattern to experience an 
Eckhaus instability; and (b) the transition from the fca pattern to the fci pattern, when the 
value of the control parameter is changed from g = 20 to g — 19, causing the fcs pattern to 
go through a saddle-node bifurcation. In both cases the results are obtained by a numerical 
integration of the seven truncated mode equations ( |11.57| l, for modes 6_3 to 63, using the 
same parameters as in Fig. |11.3| Details of the transitions are discussed in the text. From 
Kenig et al. (2009a). Copyright (2009) American Physical Society. 
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increase of g that crosses the Eckhaus instabihty hne it is simply the mode whose 
hnear growth rate is greatest that is selected. For a slow temporal ramp of the control 
parameter g = ar, with a ^ 1, they encountered a more interesting competition 
between the different patterns. They showed that as the control parameter is ramped 
a sequence of patterns start to grow one by one, yet the growth rates increase with 
each pattern that emerges. This resembles a balanced race in which the slow runners 



are allowed to start running before the fast ones. Nevertheless, Kenig et al. (2009a 



were able to predict the winning pattern, and its dependence on the ramp rate a. 
In all cases that were checked, simulations of the original equations of motion of the 
resonators (11.35) confirm the results based on the BCL amplitude equation. 



11.6 Continuous Amplitude Equations: 
Example IV - Intrinsic localized modes (ILMs) 

11.6.1 Derivation of the PDNLS equation 

As our final example we focus on a different type of nonlinear states, namely, intrinsic 
localized modes fILMs), also known as discrete breathe r s or lattice solitons fOvchiu 



nikov and Erikhman, 1982 



Sievers and Takeno, 1988 



Campbeh et al, 2004; Ma- 



niadis and Flach, 2006). These localized states are intrinsic in the sense that they 



arise from the inherent nonlinearity of the resonators, rather than from extrinsically- 
imposed disorder as in the case of Anderson localization. ILMs were observed by Sato 
et al. (|2003a[ |20036[ |2004[ |2006[ |2007[ |2008[) in driven arrays of micromechanical res- 



onators. They were also observed in a wide range of other physical systems including 



coupled arrays of Josephson junctions (Trias et al., 2000 Binder et al., 2000), coupled 



optical waveguides (Eisenberg et al., 1998 Eisenberg et al., 2001 Cheskis et al., 2003), 



two-dimensional nonlinear photonic crystals (Fleischer et al., 2003), highly- nonlinear 



atomic lattices (Swanson et al, 1999), and antiferromagnets (Schwarz et al., 1999 Sato 



and Sievers, 2004). Thus, the ability to perform a quantitative comparison between 



our theory and future experiments with large arrays of MEMS and NEMS resonators, 
may have consequences far beyond the framework of mechanical systems considered 
here. 

We follow the work of Kenig et al. (2009 6|, whose goal was to predict the actual 
physical parameters, in realistic arrays of MEMS and NEMS resonators, for which 
ILMs can form and sustain themselves. Such predictions may have practical conse- 
quences for actual applications exploiting self-localization to focus energy, and others 
that may want to avoid energy focusing, for example in cases where very large os- 
cillation amplitudes may lead to mechanical failure. Again we wish to formulate our 
analysis in terms of a continuous amplitude equation and to display the range of stable 
ILMs on a reduced diagram — as we did for extended modes in Fig. |11.3| — helping to 
describe the general qualitative behavior as physical parameters are varied. 



We start with the same form of the equations of motion ( 11.35 ) that we used in the 



preceding section, with two differences: (1) We wish to keep an explicit parameter with 
which we can vary the linear damping, thus we define the small expansion parameter 
as = £7, with e <C 1, and 7 of order unity; and (2) We use a negative sign before 
the coupling coefficient D to model elastic coupling between adjacent beams, which 
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is stronger as the separation between neighbors increases, thus acting to stiffen the 
resonators. This leads to a dispersion curve that has a positive slope, or a positive 
group velocity. The coupling mechanism in the experimental setups in which ILMs 
were observed by Sato et al. (|2003al 1200361 120041 12006| |2007| 120081) is of this kind. Our 



equations of motion then become 

Uji + ^T^n + [1 — eh cos 2ujpt ) u„ — i_D(u„+i — 2 



-ul + wlun = 0, (11.58) 



with hats to be removed later by additional scaling. 

An experimental protocol for producing ILMs in an array of resonators with a 
stiffening nonlinearity — albeit not the one we use below — is to drive the array at the 
highest-frequency extended mode. As the resonators are collectively oscillating at this 
mode, the frequency is raised further, which through the stiffening Duffing nonlinearity 
results in an increase of the oscillation amplitude up to a point at which the extended 
pattern breaks into locahzed modes (Sato et al, 2003a Sato et al., 2006). With this in 
mind — and concentrating on the case of elastic coupling where the highest-frequency 
mode uj ~ ^/T+~2D is the staggered mode, in which adjacent resonators oscillate out 
of phase — we write the displacement of the n*^ resonator as 



£i/2[^(X„,r)e^("*-™) -f c.c] -)-e3/2y(i)(t,f,l„) -f 



(11.59) 



with slow temporal and spatial variables T — et and Xn — e^^'^n. As usual, we take the 
parametric drive frequency to be close to twice uj by setting ujp = uj + eil/2, introduce 
a continuous spatial vari able X in place of A"„, and substitute the ansatz (11.59) into 
the equations of motion (11.58) term by term. Up to order e^^^ we have 



,h 




n = e'^/^iuJ^e'^^'-^''^ + C.C., 



i{ujt—7rn) 



+ 0(e*^"*,e'^'^") + c.c.. 



.3/2 (1) 
e "nil' 



(11.60a) 
(11.60b) 

(11.60c) 

(11.60d) 
(ll.eOe) 



and 



ulii^ = £3/2 .^|^|2^gz(<^i-.r«) ^ o{e'^'^\ e'3^") + C.C., 



(11.60f) 



) are fast oscillating terms with temporal frequency 3uj or spatial 



where 0{e'^'^\e'^''" 
wavenumber Stt. 

At order e-^^^ the equations of motion (11.58) are satisfied trivially. However, 
once again at order e^^^ we encounter secular terms — in this case, proportional to 
gi(wt-7rn) — g^^^ must apply a solvability condition, requiring all such terms to vanish. 
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Again, it is this condition that leads to a partial differential equation (PDE) describing 
the slow dynamics of the amplitudes of the resonators, 



2iijj 



di/j 
df 



(3 + ia;ry)|V'|V 



2 dx^ 



0. 



(11.61) 



Note that while e'('^*+^") — g«('^*-'r")^ jf -^e were to consider an arbitrary mode of 
wave number q instead of tt, the parametric term would have forced us to apply 
another solvability condition, requiring terms proportional to e*("*+5") to vanish. This 
was exactly the situation in the preceding section [11. 5[ where we were forced first to 
consider an ansatz based on counter propagating waves as the 0(e^/^) solution for w„, 
leading to a system of two coupled amplitude equations ( 11.39[ ), after which a second 
scaling was used to obtain a single amplitude equation (11.47). Here we can get away 
with a single step. 

By means of rescaling, 



(11.62) 



we transform eqn (11.61 1 into a normalized form 
. dtp 9^-0 



'ar dX^ 



ijijj - (2 + M/)|-0|^V + h-)Jj*e 



*„2iT 



(11.63) 



We then perform one final transformation ip — > tpe^'^ and arrive at an autonomous 
PDE, which is the final form of our amplitude equation, 



.dtp 



dX^ 



+ (1 - t"f)ij - (2 + ^^7)|?MV + h^j* 



(11.64) 



Equation (11.63) with 77 = is called the parametrically driven damped nonlinear 
Schrodinger equation (PDNLS). It models parametrically driven media in hydrody- 



namics ( Zhang and Vihals, 1995| |Wang and Wei, 1997 


Wang and Wei, 1998 Miao 


and Wei, 1999 


) and optics (|Longhi, 1996 


Sanchez-Morcillo et al, 2000 


), and was 



also used as an amplitude equation to study localized structures in arrays of cou- 



pled pendulums (Denardo et al, 1992 Chen, 1994 Alexeeva et al., 2000). Recently, a 



pair of linearly-coupled PDNLS equations was used to model coupled dual-core wave 



guides (Dror and Malomed, 2009). Equation (11.64) has the form of a forced complex 



Ginzburg-Landau equation (Burke et al., 2008) but with specific coefficients that are 



derived, via the scaling performed in (11.59) and (11.62), from the underlying equa- 



tions of motion (11.58). We note that considering the equations of motion (11.3) (yet 



still with a negative sign before D) as our stating point instead of eqns ( 11.58 ) leads to 



the same eqn (11.61) as above, but with slightly different coefficients. Thus, applying 
modified scahng (11.62) yields exactly the same amphtude equation (11.64). 
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11.6.2 Analyzing and solving the equation 



A remarkable feature of the amplitude equation (11.64) is that for 77 = it has exact 
steady-state solitonic solutions, as shown by Barashenkov et al. (1991), 

^±{X) = A±e-*®±sech [A± {X - Aq)] , 



where Aq is an arbitrary position of the soliton. 



1 ± ^K^ -72, and cos(2e±) = ±\ 1 



(11.65) 



(11.66) 



This pair of solitonic solutions exists for 7 < /i, above the dotted line in Fig. 11. 5| It 
was shown by Barashenkov et al. (1991) that the soliton is unstable for all values 



of 7 and h, while the ^'_|_ soliton is stable in a certain parameter range. A simple linear 
stability analysis shows that the zero solution il^{X) — 0, which exists for all parameter 
values, is stable only for h < y/l + 72. This inequality, indicated by a dashed line in 
Fig. 11.5 , also determines an upper stability limit for localized solutions of eqn ( 11.64), 



as they decay exponentially to zero on either side. 



We follow Kenig et al. (20096) in constructing an approximate analytical expression 
for the localized solution of the full amplitude equation (11.64), with 77 > 0, imple- 



menting a method introduced by Barashenkov et al. (2003). To this end, we consider 
a function of the same form as , 




Fig. 11.5 Stability diagram for localized solutions of the amplitude equation (11.641 in the 
h vs. 7 plane. The dotted line is the lower existence boundary for = 0, namely h — ^. The 
dash-dotted line is the approximate low boundary for rj = 0.1, given by eqn (11.741. Above 
the solid line the ^P^. solution of the PDNLS equation with 77 = is unstable with respect to 
a Hopf bifurcation. The dashed line is the line h = \/l + 7^ above which the zero solution is 
unstable. Red *s are points for which linear stability analysis shows that perturbations away 
from the soliton solution ^{X) grow exponentially for rj — 0.1, hence the soliton solution is 
unstable. Blue dots represent points for which the solution ^(X) is stable according to the 
linear analysis. From Kenig et al. (20096). Copyright (2009) American Physical Society. 
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ip{X, T) = a(r)e-'^('^)sech [a{T) {X - Xo)] , 



(11.67) 



except that o and 9 are now time-dependent. We multiply Eq. (11.641 by i/'*, subtract 
the complex conjugate of the resulting equation and get 



' dT dX\dX^ ^ dX 



(ii.e 



By substituting -0 — ji/ije integrating over X' = X — Xq, and assuming that ip ~^ 
and dip/dX — >■ as \X\ — > oo, we obtain a spatially- independent integral equation 



d 

dT 



J \i:\^dX' = 2 J \i;\^[hsm{2x)-l]dX' -2r] J \i:\^dX'. 



(11.69) 



Substituting the ansatz (11.67 1 into cqn ( 11.69 ), we obtain the time evolution equation 
for a 



da 
dT 



2a{hsm{2B)--f-fia^), (11.70) 

where f) = 2r]jZ. The time evolution equation for d is derived in a similar way by 
multiplying eqn ( 11.64 ) by -0* , adding the complex conjugate of the resulting equation, 
substituting the ansatz (11.67), and integrating over space to yield 

dQ 

— = /icos(26l) + 1 - 



(11.71) 



Equations (11.70) and (11.711 have the same form as the equations obtained by 



Barashenkov et al. (2003 ) , whose fixed points are 



which has to be positive, and 



1 + ry^ 



hcos{20±) = a| - 1, 
ft,sin(26'±) = 7 + rja'].. 



(11.72) 



(11.73) 

A linear analysis of these stationary points shows that (o-|-,0+) and {a^,9^) are a 
stable node and a saddle, respectively (Barashenkov et al., 2003). The saddle-node 
bifurcation point of these solutions occurs at 

1 + f] 



hsniv) 



^2 - where V ^ -ij, 



(11.74) 



as long as 7?7 < 1. This is the approximate minimal driving strength required to 
support a localized structure in the array, in the presence of linear and nonlinear 
dissipation. It is indicated by a dash-dotted line in Fig |11.5[ for 77 = 0.1. 
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The approximate stable localized solution of the amplitude equation (11.64) is 
therefore given by 

^app(^) = a+e-'«+sech(a+(X - Xo)). (11.75) 



Substituting this expression into eqn (11.59) yields an approximate expression for the 



displacements of the actual resonators in the array of the form 

/ l2eujn 



Unit) 



2\l — - — a+sech 



D 



-n - Xo 



cos (ujpt ~ nn — 6+) . (11.76) 



To obtain accurate solutions one has no choice but to solve the amplitude equation, 
or the underlying discrete equations of motion, numerically. We do so with the equa- 



tions of motion (11.58) by initiating them with the approximate expression (11.76) at 



a value of h just above the saddle node hsn (11.74). We then perform a quasistatic 



upward sweep of h, raising h in small increments and waiting for transients to de- 
cay at each step. To obtain the stationary solution of the amplitude equation we set 
dip/dT = in eqn (11.64) and solve it numerically as a boundary value problem over 
an interval of length L, with boundary conditions = 0) = i^i^ = L) — 0. We use 
the approximate expression tjjg,pp{X) [eqn (11.75 )] as an initial guess. Having identified 
an upper stability boundary h — -\/l + 7^ and an approximate lower existence bound- 



ary, given by Eq. (11.74), we must make use of the numerically obtained localized 



solutions ip{X) in order to examine the stability within these boundaries [see Kenig 
et al. (2009 &) for details]. The stability diagram of both the analytical solution ^+ 
for 77 = ( [Barash enkov et ai, 1991) and the numerical solution il^{X) for 77 = 0.1 
are displayed in Fig. |11.5| Numerical integration of the underlying equations of mo- 



tion (11.58) confirms the stability analysis, based on the amplitude equation (Kenig 



et al., 2009 &) 



Figure 11.5 highlights the effects of nonlinear damping on localized solutions. The 
first effect is to raise the lower existence boundary. This is explained by the fact that 
the additional energy lost through nonlinear damping has to be compensated by an 
increase in the strength of the parametric drive, as predicted by the approximate 



expression (11.74). The second effect is that nonlinear damping increases the area 



in the (/i, 7) parameter space where solitons are stable (blue dots). In particular, 
the shape of the unstable region for 77 > (red *s) becomes qualitatively different. 
There are values of 7 for which an increase in the drive amplitude h initially induces 
an instability of the soliton, while upon further increase of h the soliton regains its 
stability. This can be explained by noting that the amplitude of the soliton — given 



approximately by Eq. (11.72) — increases as h becomes larger, thereby enhancing the 



effect of nonlinear damping. This increase of damping exerts a similar stabilizing effect 
as that of increasing 7 in the absence of nonlinear damping. 

11.6.3 Brief survey of applications 



Kenig et al. (20096) studies a host of dynamical phenomena using their version of 



the PDNLS equation (11.64). These include questions concerning the interaction of 



pairs of solitons, which can be either attractive or repulsive, depending on the relative 
phase of the two solitons; the possibility of pairs of solitons to form bound states; and 
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the ability of a boosted soliton spontaneously to split into two. We urge the reader 
to consult the original work for additional detail, and only emphasize the exquisite 
agreement between the predictions made with the amplitude equation (11.64), and 
numerical simulations of the underlying equations of motion (11.58). 

There is one particular question that we wish to address here which deals with 
the procedure for generating solitons in actual experiments, where one cannot simply 
introduce an approximate soliton as an initial condition like one does in a simulation. 
Moreover, it is not obvious how dynamically to form solitons starting with a motionless 
array of resonators, as one needs to take the system sufficiently far from the basin of 
attraction of the zero solution ip{X) = 0, which is also stable whenever solitons are 
stable. The most direct procedure for avoiding the zero solution, starting from weak 
random noise, is to drive the system with h > -^1 + 7^, so neither the zero solution nor 
the soliton solutions are stable. As a consequence, a non-zero pattern develops. Stable 
solitons can then be formed by lowering the drive amplitude to a value h < -^/l + 7^ 
for which the zero solution and the soliton solutions are both stable, if the non-zero 
pattern that was obtained is outside the basin of attraction of the zero solution. 

This simple procedure, sometimes called self trapping — which could be imple- 
mented experimentally in a straightforward manner — is demonstrated in the top panels 
of Fig. 11.6 showing a numerical simulation of the equations of motion (11.58) with 
fixed boundary conditions, using N = 199 resonators. One can see that the initial 
transient that forms becomes unstable upon lowering the drive amplitude, giving rise 
to the formation of a number of solitons. Note that before reaching steady state, a pair 
of solitons merges into one, and another pair attracts and forms a bound state. Both 
of these effects were studied by Kenig et al. (20096). The emerging isolated solitons 
agree well with the approximate analytical form (11.76), determined earlier, with only 
their central positions Xq used as fitting parameters. 

A more controlled procedure for generating solitons would be to initiate the array in 
a particular non-zero state and then to drive it outside its known stability boundaries. 
This has been considered in the past in systems without nonlinear damping, using 
the non-zero uniform solution of the PDNLS (Barashenkov et ai, 20031. However, 
it is known for systems with 77 = that the uniform solution is always unstable 
against weak modulations and so may be difficult to access dynamically. What |Kenig| 



et al. (20096) discovered was that nonlinear damping can act to stabilize the non-zero 



spatially-uniform solution, making the procedure for generating solitons through a 
modulation instability of a uniform state possibly relevant for experiments. 

We demonstrate the use of the stable uniform solution in the dynamical forma- 
tion of solitons in the bottom panels of Fig. |11.6[ showing a numerical simulation of 
the equations of motion (11.581 with periodic boundary conditions, using N = 200 
resonators. The array is initiated with the large-amplitude uniform solution and is 
driven within the stability boundary of this state, as calculated from the amplitude 
equation (11.64). After a long time during which the uniform solution indeed remains 
stable, the drive amplitude is lowered below the stability threshold for this solution, 
but within the stability boundaries of the soliton solutions, and solitons are formed 
via a modulation of the unstable uniform state. 
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Fig. 11.6 (Color) Numerical simulation of the equations of motion (11.581 showing the 
dynamical creation of solitons, for y — 1, r; = 0.3, D — 0.25, e — 0.01, and cjp = 1.002tj. 
Plotted are the absolute values of the displacements of the resonators, which alternate between 
positive and negative values. Left panels show the time evolution, with m counting the number 
of drive periods. Right panels show the initial (black dots) and final (blue circles) states along 
with the analytical form of the solitons (green solid line), using only their central positions 
Xo as fitting parameters. Top panels: A simulation of 199 resonators with fixed boundary 
conditions is initiated with random noise and a drive amplitude of /i = 5, which is above 
the upper stability limit, h = ^/l + 7^ ~ V^, for both the zero-state and the solitons. At 
time m = 600 drive periods, after a non-zero transient (black -|-s in the right panel) has 
developed, the drive amplitude is lowered to h — 1.35 < -\/2, yielding stable solitons. Bottom 
panels: A simulation of 200 resonators with periodic boundary conditions is initiated with 
the uniform non-zero solution and a drive amplitude of h = 1.3, which is above the calculated 
stability threshold, hth — 1.26, for this state (Kenig et al. 20096). After m — 10000 drive 
periods during which the uniform state remains stable, the drive amplitude is lowered to 
h = 1.2 < hth, yielding stable solitons. From Kenig et al. (20096). Copyright (2009) American 
Physical Society. 
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